Monte Carlo study of two-dimensional Bose-Hubbard model 
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One of the most promising applications of ultracold gases in optical lattices is the possibility to 
use them as quantum emulators of more complex condensed matter systems. We provide benchmark 
calculations, based on exact quantum Monte Carlo simulations, for the emulator to be tested against. 
We report results for the ground state phase diagram of the two-dimensional Bose-Hubbard model 
at unity filling factor. We precisely trace out the critical behavior of the system and resolve the 
region of small insulating gaps, A <C J. The critical point is found to be (J/U) c = 0.05974(3), in 
perfect agreement with the high-order strong-coupling expansion method of Ref. In addition, 
we present data for the effective mass of particle and hole excitations inside the insulating phase 
and obtain the critical temperature for the superfluid-normal transition at unity filling factor. 
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In the last few years, manipulation of quantum 
gases in optical lattices has been characterized by fast 
and striking advances in trapping techniques (see e.g. 
Q and [H for a review), with the main remaining chal- 
lenge being the addressability of single sites. It is to- 
wards this direction that some experimental groups are 
devoting their efforts nowadays [5|. Access to a single 
site would enable in situ measurements of observables 
of interest and direct measurement of correlation func- 
tions, the knowledge of which would ease the study and 
characterization of new exotic states of matter. While 
ultracold gases in optical lattices are of interest on their 
own, one could also think of a broader and more ambi- 
tious project of using such systems as quantum simula- 
tors of difficult-to-solve condensed matter systems and 
models. One prominent example is quantum magnetism 
in electronic systems which may be relevant to high T c 
superconductivity. Since such systems are theoretically 
hard to address, one could alternatively think of mim- 
icking models of interest (Hubbard models, for example) 
with ultracold gases in optical lattices. 

It is in this framework that the DARPA agency has 
developed and funded a program whose goal is building 
fermionic and bosonic optical lattice emulators. There 
are two main challenges to meet: addressability of sin- 
gle sites and engineering of exchange interaction among 
atoms. The addressability of single lattice site is cru- 
cial not only for local measurements but also for ma- 
nipulation of single atoms, which would open up the 
way to applications in quantum computing. Engineering 
spin exchange interactions is essential in order to study 
quantum magnetism. It has already been shown that 
two-component boson systems with properly tailored ex- 
change interactions, can be used to realize quantum spin 
Hamiltonians @, 0, 0. Altogether, the optical lattice 
emulator, the first example of the special purpose quan- 
tum simulator, would enable one to explore new exotic 
states and answer open questions in the fields of quan- 
tum magnetism and superconductivity, including the in- 
terplay between the two (e.g. by determining ground 
states of Hamiltonians with competing orders) . 



Within the quantum gas microscope implementa- 
tion individual atoms are magnetically transported 
from a 2D surface trap in the focal plane of an ultra- 
high aperture objective to a spatially separated vacuum 
chamber Q . The most natural first step for understand- 
ing advantages and limitations of this technique of atom 
imaging is to calibrate it against the simplest correlated 
2D system, the Bose-Hubbard Hamiltonian on the square 
lattice: 
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where b\ and bi are the bosonic creation and annihilation 
operators on the site i, J is the hopping matrix element, 
U is the on-site repulsion and //, = fj, — V(i) is the dif- 
ference between the global chemical potential [i and the 
confining potential V(i). At zero temperature and inte- 
ger filling factor, the system features the superfluid(SF)- 
to-Mott-insulator(MI) phase transition with the MI 
phase being uniquely characterized by the energy gap A 
to create a particle-hole excitation. The ground state 
phase diagram of the homogeneous system (in the fi/U 
vs J/U plane) has a characteristic lobe shape with the 
system being in the MI state inside the lobe and SF out- 
side. 

In experiments, gases in optical lattices are confined 
by an external potential. So far, this has resulted in lim- 
itations in the observation of a quantum phase transi- 
tion due to measurement averaging over the whole cloud. 
With the high-resolution quantum gas microscope, mea- 
surements can be performed locally and averaging over 
the inhomogeneous system can be avoided. The first 
goal of the bosonic quantum emulator is to map out the 
ground state phase diagram. The standard approach is 
based on a local chemical potential approximation where 
the density of the homogeneous system with the chemical 
potential 



^ = „ - V{i) , 



(2) 



is identified with the density at the site i of the inhomo- 
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FIG. 1: (Color online) Phase diagram of the first MI-SF lobe. 
Solid circles are numerical data, with error bars shown but 
barely visible. The inset is a blow up of the region close to 
the tip. Dashed lines represent the critical region as calculated 
from finite size scaling. 



geneous system. 

In this paper we provide benchmarks for the bosonic 
quantum emulator to be tested against. We report re- 
sults of large-scale exact quantum Monte Carlo simula- 
tions of model [TJ by the Worm algorithm [TlJ . We focus 
on the homogeneous case and unity filling factor. Worm 
algorithm allows efficient sampling of the single-particle 
Green function. Precise data for the Green function en- 
able us to carefully trace out the critical behavior of the 
system and resolve the phase diagram in the region of 
small insulating gaps, A <C J. We also present data for 
the effective mass of particle and hole excitations inside 
the insulating phase. Effective masses characterize the 
phase transition away from the tip of the lobe. Here the 
transition is described by the physics of the weakly in- 
teracting Bose gas in the limit of vanishing density 

In order to completely characterize the sys- 
tem the full phase diagram in the parameter space 
(fi/U, J/U, T/J), where T is the temperature, is needed. 
Here we limit ourselves to studying ground state proper- 
ties and calculating the critical temperature for the SF- 
normal transition at unity filling factor. An exhaustive 
finite temperature study of the system is in progress in 
another group [T^. 

We now turn to the presentation of our results. 
The procedure used to determine the ground state phase 
diagram and extract effective masses of particle and hole 
excitations from the Green function was discussed in de- 
tails in Ref. [HI]. In Fig. [1] we present results for the 
ground state phase diagram corresponding to unity fill- 
ing. The inset shows the region around the tip. Circles 
represent the simulation data while dash lines are ob- 
tained from the finite size scaling analysis. Simulations 
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FIG. 2: (Color online) Finite size scaling of the energy gap at 
the tip of the lobe. Lines represent linear fits used to extract 
the critical point. The critical point can be directly read from 
the intersection of the curves: (J/U) c = 0.05974(3). 

were done for linear system sizes L — 10, 20, 40, 80. We 
do not see any significant size effect up to J/U ~ 0.057. 
In order to extract the position of the critical point at 
the tip of the lobe and determine the extension of the 
critical region, the standard finite size scaling argument 
was used (see Ref. 13]), with the critical exponent for 
the correlation length v = 0.6715. The finite size scal- 
ing of the energy gap is presented in Fig. [2] One can 
directly read the position of the critical point from the 
intersection of the curves: 

(J/U) c = 0.05974(3) (n = I). (3) 

Equation ([3]) and Fig. [T] constitute the most precise 
quantum Monte Carlo simulation for the Hamiltonian [T] 
which is in perfect agreement with the result of Ref. [f[ , 
where the authors carried out a strong coupling expan- 
sion up to 13-th order. Note that the critical region in 
Fig.[T]is resolved with accuracy <C J, i.e. for gaps A < J, 
which is crucial for studies of the emerging relativistic 
physics at the lobe tip. 

In Fig. [3] we plot effective masses for particle (cir- 
cles) and hole (squares) excitations. Dispersion rela- 
tions were fitted by a parabola, with the exception for 
J/U — 0.059 where we used a relativistic dispersion re- 
lation. Close to the tip of the diagram, the action is 
isotropic in space and imaginary time, giving rise to a rel- 
ativistic behavior [io| . In the limit J/U — ► 0, where one 
can calculate effective masses perturbatively, our data 
converge to the analytical result (dashed lines). To the 
first order, the perturbative expansions are given by (we 
set the lattice pariod and Planck's constant equal to 
unity): 

Jm + = 0.25 - 2 J/U, Jm^ = 0.5 - 8 J/U (4) 
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FIG. 3: (Color online) Effective mass for particle (circles) 
and hole (squares) excitations as a function of J/U. The 
exact results at J/U = are m+ = 0.25/ J and m_ = 0.5/ J. 
By dashed lines we show the lowest order in J/U correction 
to the effective masses. Close to the critical point the two 
curves overlap, directly demonstrating the emergence of the 
particle-hole symmetry. At J/U — 0.059, the sound velocity 
is c/J = 4.8 ±0.2. 

for particle and hole excitations respectively. On ap- 
proach to the critical point, instead, data for particle 
and hole excitations are merging together, in agreement 
with the emergent particle-hole symmetry in the criti- 
cal region. From the fit done at J/U — 0.059 using the 
relativistic dispersion relation, we extract the value of 
the sound velocity c/J = 4.8 ± 0.2 and effective mass 
to* = 0.015 ± 0.0015. We would like to point out that 
effective masses can also be extracted using the method 
of Ref. [l[, although we did not find any calculation in 
the literature. 

In Fig. |4] we show the phase diagram for the SF- 
normal transition at integer filling factor n = 1. The 
transition is of Berezinskii-Kosterlitz-Thouless type [lij . 
The critical temperatures were found from extrapolation 
to infinite system size of the standard finite size scaling 
for the Kosterlitz-Thouless transition (see e.g. Ref. [l5j). 
In the figure, circles are numerical results and dashed 
lines are analytical expressions in the two limiting cases. 
In the weakly interacting regime the critical temperature 
is given by: 

T (5) 
to In^/mU) ' 

where n is the density (n = 1 in this case), to is the 
mass (m = 1/2 J in the lattice), and £ is a dimension- 
less parameter which was found numerically in Ref. (l5| 
to be £ = 380 ± 3. On the approach of the critical 
point, instead, one can use the following scaling argu- 
ment. Close to the critical point the superfiuid density is 
p s (T = 0) ~ C" 1 ~ t", where t = (U/J-(U/J) C ). Under 
the assumption p s (T — 0) ~ p{T c ), which should hold for 



FIG. 4: (Color online) Finite-temperature phase diagram at 
filling factor n = 1. Solid circles are simulation results (the 
dotted line is to guide an eye), error bars are plotted. Dashed 
lines are analytical results for the weakly interacting gas [see 
Eq. {5}] and for the strongly interacting gas close (J/U) c . 



low enough critical temperatures, i.e. close to the critical 
point (U/J)c, one concludes that T c ~ p(T c ) ~ t v . The 
dashed line in the plot is a fit done using the function 
fix] = At" , where A=0.49(2) is a fitting parameter. On 
both sides, numerical results clearly converge to the an- 
alytical expressions. 

Many interesting phenomena happening at zero or 
nearly zero temperature have not been observed yet. This 
is because, so far, it has been a challenge for experimen- 
talists to reach low enough temperatures. In order to 
overcome this challenge, one can exploit the inhomogene- 
ity of the entropy distribution of the harmonically and 
optically trapped gas. The idea has been originally pro- 
posed in Ref. [r| and [13], where the authors suggest 
several cooling protocols. Some of the protocols make 
use of the filtering scheme of Ref. [l8[ , others require spin 
dependent lattices. In all the protocols the main idea is 
to relocate the entropy by removing a small fraction of 
the particles carrying almost all the entropy. Recently 
it has also been suggested a cooling protocol based on 
coupling entropic particles with a system at lower tem- 
perature (i.e. a "refridgerator" ) [lj|. Having in mind a 
setup, similar to the one described in Ref. [l(| and [l7l ]. 
we would like to suggest a simple and efficient cooling 
protocol which does not require coupling to a refridger- 
ator or exciting particles to a different internal energy 
level. Consider a system with U 3> |/ij D — p,i=o\,J,T 
and — 3> J, where ^ is given by Eq. ^ty, 

i = is at the potential minimum while i ~ io corre- 
sponds to the boundary of the system. The condition 
U \jii — /ii=o\,J means that the groundstate density 
is essentially uniform. The condition U 3> T guaranties 
that the entropy is located in a thin peripheral region of 
the system. In view of the condition |^ — 3> J, the 
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FIG. 5: Sketch of the proposed cooling protocol. On top, a 
shallow and a deep trap are superimposed. At the bottom, 
the corresponding densities of states are sketched. The dashed 
line represents the Fermi level. Adiabatically displacing the 
shallow trap results in displacing entropic particles with con- 
sequent dramatic decrease of entropy per particle inside the 
deep trap. 

elementary excitations of the system are single-site (lo- 
calized) particles and holes obeying Fermi statistics (in 
the real space) [l6j . 

Take two superimposed traps, a very steep one and 
a very shallow one. The density of energy levels of the 
shallow trap is very high as compared to that of the steep 
one. The shallow trap thus carries most of the entropy. 
Adiabatically displacing the shallow trap would then re- 



sult in displacing "entropic" particles and ultimately sep- 
arating them from the particles of the steep trap. Aca- 
demically speaking, the very last stage of this process 
cannot be adiabatic in view of the exponentially sup- 
pressed tunnelling between two traps. However, this non- 
adiabaticity simply means that the two systems become 
independent, so that the (now dramatically decreased) 
entropy of the steep trap has nothing to do any longer 
with the state of the particles in the shallow trap. The 
procedure is sketched in Fig. [51 For better results, one 
can perform a number of cooling cycles consisting of the 
above-described entropy relocation procedure followed by 
adiabatically increasing (and subsequent decreasing) the 
J/U ratio in order to lift the localization constraint pre- 
venting redistribution of tiny fraction of remanent parti- 
cle and hole excitations in the bulk towards the perime- 
ter. Another possibility for entropy relocation consists of 
progressively lowering the confinement below the Fermi 
level in order to evaporate the entropic particles from the 
system. 

In conclusion, we have presented numerical results 
for the single component Bose-Hubbard model. We have 
confirmed previous calculations 1] for the critical point 
at the tip of the n = 1 lobe and presented results for 
the effective masses for particle and hole excitations. 
We have determined the critical temperature for the SF- 
normal transition for the unity filling case. These bench- 
mark calculations provide a robust test for the bosonic 
optical lattice emulator. 

We are grateful to Tin-Lun Ho for useful discus- 
sions. This work has been supported by the DARPA 
OLE program. 
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